Efficiency analysis of reaction rate calculation methods using analytical models I: The 

2D sharp barrier 
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We analyze the efficiency of different methods for the calculation of reaction rates in the case of a 
simple 2D analytical benchmark system. Two classes of methods are considered: the first are based 
on the free energy calculation along a reaction coordinate and the calculation of the transmission 
coefficient, the second on the sampling of dynamical pathways. We give scaling rules for how this 
efficiency depends on barrier height and width, and we hand out simple optimization rules for the 
method-specific parameters. We show that the path sampling methods, using the transition interface 
sampling technique, become exceedingly more efficient than the others when the reaction coordinate 
is not the optimal one. 

PACS numbers: 



I. INTRODUCTION 



As Molecular dynamics (MD) is limited to microscopic 
systems and time scales, most chemical or biological re- 
actions can not be simulated using straightforward MD. 
One can literally wait ages before detecting a single event 
in a typical computer simulation. In the early 1930s, 
Wigner and Eyring made the first attempts to overcome 
this problem by introducing the concept of the Transi- 
tion state (TS) and the so-called TS Theory (TST) ap- 
proximationiiS. Later on, Keck^ demonstrated how to 
calculate the dynamical correction, the transmission co- 
efficient. This work has later been extended by Ben- 
nett, Chandler"^ and othersSii, resulting in a two-step 
approach. First the free energy as function of a reac- 
tion coordinate (RC) is determined. This can be done 
by e.g. Umbrella Sampling (US)^ or Thermodynamic In- 
tegration (TI)^. Then, the maximum of this free energy 
profile defines the approximate TS dividing surface and 
the transmission coefficient can be calculated by releasing 
dynamical trajectories from the top. This approach is, 
in principle, exact and independent of the choice of RC. 
However, the method becomes inefficient when the trans- 
mission coefficient is small. A proper choice of the RC 
can maximize the transmission coefficient and is hence 
crucial for the efficiency of the method. 

There exist different formalisms for the transmission 
coefficient formula which differ in the way trajectories 
are counted. We discuss the standard Bennett-Chandler 
{BC% the history dependent BC {BC2% and the effec- 
tive positive flux (EPF )i94l - formalism. We show that the 
latter should always be preferred due to a lower average 
pathlength and a faster convergence. However, whenever 
a lot of correlated recrossings occur, the transmission co- 
efficient will be very low and all these methods become 
inefficient. In high dimensional complex systems it can 
be a very difficult task to find a proper RC. Moreover, 
whenever the dynamics is diffusive, even an optimal RC 
can result in a very low transmission and hence a poor 
efficiency. 



A new approach came with Transition Path Sampling 
(TPS)i^ that is not based on the free energy barrier as 
starting point. TPS is rather an importance sampling 
of dynamical trajectories. Hence, it is a Monte Carlo 
(MC) sampling in path space rather than phase space. 
The TPS method has been advocated as a method that 
does not need a RC and is akin to 'throwing ropes in 
the dark—. This might be true if one wants to sample 
a set of reactive trajectories, but it is not for the calcu- 
lation of reaction rates. In fact, the original approach 
to calculate reaction rates within the framework of TPS 
required the definition of an order parameter and the cal- 
culation of the reversible work when the endpoint of the 
path is dragged along this parameter. For the sampling 
of reactive pathways, the order parameter needs only to 
distinguish between the two stable states. However, in 
the algorithmic procedure to calculate reaction rates with 
TPS, the order parameter becomes very similar to a RC. 
Still, it has been speculated that this approach is less 
sensitive to the problems related to an improper RC (or 
order parameter). Indeed, in this article we prove for the 
first time that this is true using the approach of Transi- 
tion Interface Sampling (TIS)i^. TIS increases the effi- 
ciency of the original TPS rate calculation considerably 
by allowing the pathlength to vary and by counting only 
positive effective crossings. The overall reaction rate in 
TIS is obtained from an importance sampling technique 
that uses a discrete set of interfaces between the sta- 
ble states. Hence, TIS could be considered a dynamical 
analogue of US in path space. For diffusive systems the 
partial path TIS (PPTIS) was invented that uses the as- 
sumption of memory lossi^. In this article we discuss the 
case of sharp barriers. Here recrossings occur mainly due 
to the wrong choice of RC. In a follow-up article we will 
treat the diffusive case. 

Up to now, it is not clear how these methods compare 
in efficiency and the need for benchmark systems has 
been put forward several times. It is not always easy to 
perform comparative calculations since it is not simple to 
know if each method is equally optimized for a specific 
system. Therefore, in this paper we analyze a system 
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for which the efficiency of the methods can be calculated 
analytically with only a few approximations. This does 
not only give a transparent comparison of the efficiency of 
the different methods, but also allows to obtain scaling 
laws for how this efficiency changes as function of the 
barrier height and width. Moreover, we give some rules 
for how the methods can be optimized, for instance, by 
choosing the proper width of the US windows and the 
position of interfaces in TIS. These rules are important 
for the simulation community, as they can be used as a 
rule of thumb in daily practice, when any method needs 
to be optimized. 

The principal component to measure the efficiency of 
the methods will be the CPU efficiency time ToS which is 
the lowest computational cost needed to obtain an overall 
statistical error equal to one. We give a detailed analysis 
of how Toff can be calculated for some very general cases 
in the appendix sections ^ and ^ It also gives the im- 
portant result of how one should divide a total fixed sim- 
ulation time over a set of different simulations to obtain 
the best overall efficiency. In Sec. we outline the first 
class of methods, the reactive flux (RF) methods, and 
present their principal formulas. In Sec. IIIII we do the 
same for the second class of methods, the path sampling 
methods. Sec. IIVI introduces the 2D benchmark system 
where the angle indicates how far the chosen RC is de- 
viated from the optimal one. Sec. is the main section 
of our paper in which we apply the different methods 
of Sec. in] and IIIII to the analytical benchmark system 
of Sec. IIVI Finally, we discuss the important point of 
hysteresis for the two types of methods and show that 
this is less likely to occur for the path sampling methods. 
We summarize the results in Sec. lVIII Moreover, to sup- 
port the readability of this paper we have added a list of 
symbols and abbreviations in App.lElandlFl 



II. FIRST CLASS OF METHODS: REACTIVE 
FLUX METHODS 

A. General formalism for RF methods 

In all combined free energy and transmission coefficient 
based methods, the rate equation follows from k ~ k{t') 
where the reaction rate k is expressed as a quasi-plateau 
value at a time t' of a time dependent reaction rate func- 
tion k{t). This function is given by the corrected flux 
through an hypersurface {x|A(a;) = A*}, that is is a col- 
lection of phase points x, defined by the reaction coor- 
dinate A(a;) and transition state (TS) value A*. The TS 
value A* is standardly taken as the maximum in the free 
energy profile along A. Both TST, BC, BC2 and EPF 
can be expressed as 



k{t) 



{\{xo)6{\{x^)-\*)x[X,t]) 
{9{\*-\{x,))) 



(1) 



where 5{-) is the Dirac delta function, 0{-) the Heaviside 
step function, xt is the phase point at a time i, and X 



is a trajectory that includes xq. Ensemble averages (. . .) 
in phase and path space are defined in Add. I A ll Eq. 
measures the flux contributed by pathways leaving the 
surface A* at t = under the influence of a correction 
functional x[A', t]. The functional x has different forms 
for TST, BC, BC2, and EPF. The rate equation (jTJ can 
be rewritten as a product of two factors: the probability 
to be on top of the barrier times the transmission function 

n{t): 



with 



and 



Here, 
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k{t) ^ PA{\*)n{t) 

{S{X{xo) - \*)) 



{0{X*-X{xo)))- 



n{t) = (A(a;o)x[^,i]>5(A(2;o)-A*)- 



(2) 



(3) 



(4) 



IS(\(xo)-\*) 



• ) implies that the ensemble {xq\ is 



constrained at the surface {xo|A(a;o) = A*}. Substitution 
of Eqs. H4I3|I in Eq. (0) using Eq. (|A5p gives back Eq. iQJ. 

Eqs. (|2I3|I show the two-step procedure. The proba- 
bility Pa (A*) and the time dependent transmission func- 
tion TZ{t) are calculated in two separate simulations. As 
for the rate fc, the unnormalizcd transmission coefficient 
TZ follows from a plateau in this time dependent func- 
tion: TZ = TZ{t'). This factor corrects for the corre- 
lated recrossings. In III Bl we discuss the methods to 
compute Pa (A*) or, equivalently, the free energy bar- 
rier AP = -A:BTlnPA(A*). Then, in HTO wc discuss the 
methods to determine the transmission coefficient TZ. 



B. Free energy methods 

1. US using rectangular windows 

Define the following 2M + 2 box functions: 

w,{x) = Q\\{x) - Afl -ir]6i[Afl + ir + 7- A(x)], 
wuix) = B\\{x) +d\- X*]e[X* - X{x)], 
Wo{x) = 9[X* - Xix)l (5) 

w,{x) = e[xix) -Xr- u - i)T]e[XH + jT + 7 - x{x)], 

with < i < M and < j < M. Here A* is the TS, Xr is 
a value in the reactant well and dX is a small length scale. 
7 and r represent the dimensions of the US windows; 
r -|- 7 is the width of the window and 7 is the overlap 
such that MT + dX = X* - Xr (See Fig. [TJ). Neglecting 
higher orders terms in dX, we can write 



Pa (A*) 



1 



(wm) 



1 



dX {0{X* - A)) dA 



[Wm)Wo- 



(6) 



To calculate Eq. ©, we can simply run an MD simula- 
tion and count the number of times that the transition 



3 




FIG. 1: Illustration of US using rectangular windows for the 
case M = 3. is a surface in the reactant well, A* is the 
TS. The binary functions Wi and Wi are depicted below the 
free energy plot where the gray areas indicate where Wi , Wi 
equal unity. 



state region interval is visited. The weight function Wq 
in the ensemble acts like an infinite waU at A = A* and 
prevents the unnecessary exploration of the product re- 
gion. However, as {wm)wo vanishingly smaU for high 
barriers, this straight-forward method will usually fail. 

Using Eq. (|A5(I and the relations WgWs = Ws, 
Ws-iWg = Ws-i for all s we can rewrite Eq. © as 

The final property is now calculated from a series of sim- 
ulations in which each pair {ws)ws , {ws^i)ws ^ {wm)wo 
so that they can be determined accurately. The im- 
plementation of US using rectangular windows via MC 
is straightforward. The standard MC sampling is per- 
formed starting from a point inside the window. As 
soon as the MC procedure generates a point outside this 
window, this point is automatically rejected and the old 
point is kept. If the procedure is performed by means of 
MD, the window boundaries simply act as infinitely hard 
walls. However, due to practical problems related to a 
discontinuous force profile, MD simulations are usually 
performed with parabolic windows instead of rectangu- 
lar ones. 



that the error does not propagate as in the case of several 
windows. On the other hand, one needs to have already 
a good idea about the shape of the barrier to construct a 
good bias f2. The series of rectangular windows is a more 
robust way to explore the barrier region when no knowl- 
edge is available. The algorithmic procedure to sample 
points in the biased ensemble (. . ■)nWo is explained in 
Sec.im 

3. Other free energy methods 

Many other methods exist for the calculation of free 
energy surfaces. TI is of equally importance and is based 
on the constrained sampling at surfaces from which the 
free energy's derivative can be obtained at a given value 
of A. Integration of these derivatives results in the re- 
quested free energy profile. In addition, many variations 
of US have been devised, such as Wang-Landau sam- 
plingi^, meta dynamicsii, and floodingi^, where the op- 
timal biasing potential is created on the fly. 

C. Transmission coefflcient calculation 

1. TST approximation 

111 TST, all pathways leading towards the product site 
B are assumed to stay in B for a very long time. The 
TST approximation uses Eq. with 

x'^^[X,t]=e{X{xo)). (9) 

If the barrier is sharp and a proper RC is taken, TST is 
a very good approximation or can even be exact. The 
calculation of TZ, in TST, requires a simple numerical or 
analytical calculation. For instance, suppose that A is a 
simple Cartesian coordinate with mass m, then A(a;o) = v 
and 

^ L^^dwwe-^i™"' l//3m 1 , , 

TZ = — = — = -^=^= (10) 

J^^dve-'^i'^^' ^2Tr /f3m VW^' ^ ' 

Therefore, the free energy profile is sufficient to obtain k 
whenever TST is valid. 



2. US using a single biasing potential 



2. BC formalism 



Instead of performing several simulations using rectan- 
gular windows, one can also use a single biasing function 



Pa{X* 



1 (wM^ ^)nwo 
AX {n-^)nwo ' 



(8) 



Again, the equivalence between Eq. ® and Eq. ijSJ) is 
easy to prove via Eq. ljA5p . The advantage of Eq. (jSJ is 



The BC equation is obtained using 

X^^[X,i]=0[A(xO-A*] 



(11) 



in Eq. The evaluation of 'R,{t) in Eq. 10} consists 

of releasing many trajectories that start from the top of 
the barrier. These trajectories only make a contribution 
different from zero if they end up at the product side of 
the barrier. It is important to note that a trajectory with 
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X{xq) < that leaves the surface heading to the reactant 
state A but finaUy ends up at the product state B at 
a time t, gives a negative contribution. The final value, 
which has to be positive per definition, results from a 
cancellation of positive and negative terms. 



3. BC2 formalism 



The BC2 equation uses 

x^^2[x,t] = e[\{xt) - \*]e[x* - \{x_t) 



t=-t' 



t=0 



t=t 



(12) 



Here, besides ending in the product state, the trajecto- 
ries integrated backward in time also have to end in the 
reactant state to give a nonzero contribution to Eq. Q). 
However, here as well, the contribution of some trajecto- 
ries arc negative. This happens when the systems starts 
with a negative A(a;o), but the forward and backward tra- 
jectories end up in the product B and reactant state A, 
respectively. These so-called S-curves (trajectories that 
cross the TS surface more than 2 times within a short 
time) are less likely to occur for sharp barriers. 



4- EPF formalism 



FIG. 2: Explanation of the transmission coefEcient 
method The dot on the surface A* indicates the point 

of release ait = Q. The arrow indicates the direction of A(a;o). 
Trajectories are integrated forward up to t = t' and backward 
in time up to t = —t'. From top to bottom, A(a::o)x using BC 
Eq. (|TTJ results in {v,Q,v,v,v,Q, -v, -v), the BC2 Eq. itHl 
gives {v,0,0,v,v,0,0,—v), and the EPF expression Eq. (I13II 
gives (v, 0, 0, 0, V, 0, 0, 0) with v = |A(a;o)| the absolute velocity 
at the initial point. 



The EPF equation arises when we take x[X, <] 



X 



EPF 



[X, t] in Eq. Q with 



EPF 



X, t] 



'O, ifA(a;o)<0, 

or if At' > A* for any t' e [-t : 0], 
or if A(.Ti) < A*, 
1, otherwise. 

(13) 

Despite the apparent mathematical more complicated 
form of Eq. (fT!^ compared to Eqs. l(TT|) and ifT^ . the 
relation is quite is natural. The approach counts only 
the first crossings and only when they are in the posi- 
tive direction (See Fig. ^ . In the effective flux formal- 
ism II13I) . all contributions are either zero or positive. In 
Ref. [ll|, we presented a similar transmission coefficient 
formula. However, in this approach the pathway was 
stopped whenever the stable state regions A and B were 
reached. For a formal mathematical proof of the equiva- 
lence between BC- and EPF- type equations see Ref. [Toj . 



5. Other transmission coefficient methods 

Some other relations for the transmission different from 
Eq. H11I13|I have been proposed. For instance, Hummer— 
proposed a relation that uses all the crossing velocities, 
in case of correlated recrossings, instead of just A(a;o). 
Ruiz-Montero et alm^ devised a transition zone method 
for diffusive barrier crossings. This will be treated in 
more detail in the follow up article that will discuss the 
diffusive barrier case. 



III. SECOND CLASS OF METHODS 
SAMPLING METHODS 



PATH 



A. General formalism of TIS methods 

Suppose Xa and Xb > Xa are such that any x with 
A(x) < A^ is at the reactant side A and any x with A(a:) > 
Xb is at the product side B of the unknown optimal TS 
dividing surface. It is important to note that A(a;) does 
not have to be a proper RC to fulfill this criterion, but 
only has to distinguish between the two stable states. 
The TIS expression is now given by 



kAi 



{h 



A) 



-Va{Xb\Xa). 



(14) 



Here (j)A gives the flux through interface Xa and is a 
history dependent function describing whether the sys- 
tem was more recently in ^ or in i?: = 1, for a phase 
point x and its corresponding trajectory, if the system 
was more recently in A than in B and hj^ = other- 
wise. In practice, the whole factor -^^^ is calculated by 
counting the number of crossings in a straight-forward 
MD simulation, divided by the number of cycles with 
hA ~ 1, divided by the time step At. The calculation of 
this factor is very cheap as interface A^i is very close to 
the basin of attraction of state A. On the other hand, 
the crossing probability 7^a(A_b|Aa) is a very small num- 
ber and can not be accurately determined by straight- 
forward MD. This is the probability that whenever the 
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system crosses interface Aa, it will cross interface \b be- 
fore crossing interface again. As interface Xb is an 
interface at the other side of the barrier, this probability 
is very small. The TIS method overcomes this by defin- 
ing Al — 1 interfaces between Xa = Aq and Xb = Aa/. 
Then, the total crossing probability can be expressed in 
a formula that contains conditional short-distance cross- 
ing probabiliticsii 

M-l 

Va{Xb\Xa) = Pa(Am|Ao) - n 'PAiK+ilXs). (15) 

s=Q 

Here, Va^Xs+iIXs) is a generalization of the previously 
described overall crossing probability and denotes the 
conditional probability that, whenever the system leaves 
state A by crossing Xa and crosses As in turn, it will 
also cross A^+i before returning to A. If the dis- 
tances As+i — Xs arc sufficiently small, the probabilities 
7^a(As+i|As) will be large enough so that they can be 
computed using a shooting algorithmS^. The shooting 
algorithm takes a random time slice from the old exist- 
ing path and makes a slight randomized modification of 
this phase point (usually only the momenta are changed) . 
Then, this new phase point is used to propagate forward 
and backward in time yielding a new trajectory. In TIS, 
this propagation of the trajectory is stopped whenever 
the system enters A or B or, equivalently, whenever Aq 
or Xm are crossed. The pathway is accepted only if the 
backward trajectory ends in A and the total trajectory 
has at least one crossing with As. The fraction of these 
paths that cross As+i as well, determines 7^a(As+i |As). 
Although the form of Eq. H15|l deceptively resembles a 
Markovian factorization, no assumption has been made. 
As 7^yi(As+i |As) are not simple hopping probabilities, but 
incorporate the full history-dependence from Xa to As, 
the relation is actually non-Markovian and exac^. 



B. TIS biasing on Amax 

In analogy with US, instead of using a discrete set of 
interfaces, we could also bias the trajectory in a contin- 
uous way using a bias on Amai^. First we can write 

^a(Aa/|Ao) - (0(A„,axm - Xb)). (16) 

Here Amax = max{A(xt)|a;t G X} where the path X is 
terminated when it leaves A and enters region A or B. 
Then, as in Eq. (|HJ) we can write 

(^'(Amax[-''^] — A_B)il^-^(Ai„ax))o 



VAiXnlXo) 



17) 



The algorithmic procedure would be the same as TIS 
without the interface crossing constraint. Instead, the 
acceptance-rejection criterion is adjusted as explained in 
Sec. lA II A continuous bias has been applied within the 
original TPS scheme. In Ref. a bias on the end point 
of the path was used. Alternatively, one can also bias the 
direction of the momenta change at the shooting point 
as was done in Ref. l25l. 



C. Other path sampling methods 



The original TPS rate calculation algorithm was in- 
troduced in Refs. It first creates ensemble of 
reactive trajectories of a fixed length. These trajectories 
should constitute a time interval that is longer than f 
where k{t) reach its plateau. Then, a second series of 
simulations is required that combines the MC of path- 
ways with a US technique. Also these simulations use a 
fixed pathlength, but these can be shorter using a rescal- 
ing approach2&. In the end, the final rate constant can be 
constructed from the results of these two types of simu- 
lations. Moroni showed that this algorithm is always less 
efficient than the TIS technique and that the computa- 
tional advancement of TIS is at least a factor two^^. The 
TIS improvement is partly due to the flexible pathlength 
so that each individual path can be limited to its strictly 
necessary minimum. Moreover, TIS has abandoned the 
shifting moves and has a stronger convergence (no can- 
cellation positive and negative terms) . Depending on the 
system and the required accuracy, the TIS approach can 
easily become more than an order of magnitude faster 
than the original approach. 

The PPTIS method reduces the average pathlength 
even further as compared to TIS using the assumption 
of memory-lossi^. In this method, the trajectories do 
not have to start at Xa- Instead, an ensemble of tra- 
jectories is generated that start and end at either As+i 
or As-i, but must have at least one crossing with the 
middle interface As. From this, four crossing probabil- 
ities can be constructed that still inhibit some history 
dependence. Once these are known for each s, the fi- 
nal overall crossing probability can be constructed via a 
recursive relationi^. The Milestoning method^! is very 
similar to PPTIS, but relies on a stronger Markovian 
assumption that the system remains in an equilibrium 
distribution at each interface. On the other hand. Mile- 
stoning uses time-dependent hopping probabilities which 
supplies a very general way to coarse-grain a dynamical 
system. The two aspects of PPTIS and Milestoning could 
be combined as was suggested in Ref. p^ . 

Finally, we mention Forward Flux Sampling 
(FFS)22ik. FES uses the same rate equation lO 
as TIS, but the sampling is different. In FES^Siffi, the 
crossing points at the next interface of the 'successful 
pathways' are stored. The next simulation uses this 
set of points to initiate new pathways. The advantage 
is that FFS does not require any knowledge on the 
distribution p{x). This allows to treat non-equilibrium 
systems as well. Moreover, EES creates effectively more 
pathways than TIS with the same amount of MD steps 
and does not rely on a Markovian assumption as in 
PPTIS. However, the correlation between the different 
pathways is much stronger than in TIS or PPTIS. 
Therefore. FFS can only be applied when the process is 
sufficiently stochastic. 
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IV. ANALYTICAL 2D BENCHMARK SYSTEM 



We consider the following 2D potential: 



oo 



if I Aj; — A* I > Rx, 
Vbar(A_L) otherwise 



(18) 



with 



A±(A,, Ay) = (A, ~ a;) cos - {Xy - a;) sing (19) 



and 



V. fA)-i°' if|A|>il^, 



(20) 



Here, H is the height of our barrier. W is the barrier 
width and Rx , Ry arc the dimensions of the reactant re- 
gion (See Fig. UJ. The chosen RC is A^;, while Xy repre- 



Rx 



Rv 



FIG. 3: The two-dimensional potential given by Eq. (I18II . H 
and W are the barrier height and width. Rx and Ry are 
the dimensions of the reactant (and product) well. A^ is the 
assumed RC, while Ay is another important degree of freedom. 
AJ and Xy are the maxima of the free energy functions for 
these coordinates. 6 is the angle between the chosen reaction 
coordinate Xx and the optimal one X±{Xx, Xy) of Eq. 1191 . 



sents another important degrees of freedom. A_L(Aa;,Ay) 
is the unknown optimal RC as its direction is orthogonal 
to the barrier ridge, the exact TS dividing surface. The 
angle 6 is, hence, a measure of the quality of the chosen 
RC. 

To facilitate the analytics we assume a simplified dy- 
namics: The particles propagate without dissipation and 
move only along the A^; direction. Once they collide 
with the walls, they obtain a new random Ay coordi- 
nate £ [Ay — i?y. A* -I- Ry] and velocities A^, taken from a 
Maxwellian distribution. This type of dynamics satisfies 
ergodicity and ensures that the reaction rate is well de- 
fined for high barriers. Although the low dimensionality 
and this dynamics might seem artificial, this minimal- 
istic model already illustrates clearly the problems that 
occur in complex systems when no adequate RC can be 



found as we will see in the forth-coming sections. More- 
over, the potential can be viewed as a projection of a 
high-dimensional complex system. In that case, Eq. 118|) 
represents a free energy in which A^, and (the unknown 
but important coordinate) Xy can be any non-linear func- 
tion of all Cartesian coordinates. For instance, Aa; could 
be the simple distance between two atoms to describe the 
formation or breaking of a bond, while Ay represents a 
complex solvent rearrangement. 

The surface potential energy at the barrier A* equals 

ViXl,Xy) = ^bardAy - A;| sine) = H(1 - 2\Xy - Xl\sin9/W) 



(21) 



and, hence, 
Pa{X* 



1 



2RxRy JX'-^-Ry 

-I3H I pipHRy siiiB/W 



dAy e-'5^(^-^«) 



2(3HRxRy sine /W 



(22) 



Here, we assumed that Rx ^ W. The transmission 
coefficient can be obtained using Eq. (0} and x^*""^ = 
^EPF (|X2I13|) . The two equations are identical 

since there are no trajectories that can cross the surface 
{x|Aa;(x) = A*} more than two times before colliding 
with the outer walls. Eqs. (|12I13|I have only non-zero 
contributions whenever the backward and forward tra- 
jectory end, respectively, at the left and right side of the 
barrier after a short time f, which is deterministically 
determined by X^ixo)- Hence, xi^it'] = xiXx, Xy)\t=o = 
9[Xx - ^2{H -V{X*,Xy))/m], which gives 



JdXyU'iXy 



JdXye-PnK^y) 



(23) 



with 

n'iXy) 

Hence, 

n = 



rOO 



V)/ri 



■ dv 1 



-0{H-V{\:,\y)) 



2PH Ry sine /W ^^,0HRy 



V27r/3m 



le/w 



which yields by Eq. ^ 



1 



(24) 



(25) 



(26) 



^/2^i#^ Rx 

The reaction rate is, thus, independent of the angle 0. 

V. EFFICIENCY SCALING 

To quantify the efficiency of the different methods, we 
will calculate the CPU efficiency time, t^s, that is de- 
fined as the minimal computational cost to obtain an 
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overall relative error equal to one. The efficiency is some- 
times defined as the inverse T~g of this quant it jiSiiS. In 
App. we show how this quantity can be computed for 
some very generic cases. As the force calculations are the 
predominant steps in any ab initio or large scale classical 
simulations, all CPU values will be expressed as an in- 
teger representing the number of required MD steps. In 
this way, we obtain a measure that is independent of the 
computational resources. 

In the following, we make two assumptions: 

• The correlation number is assumed to be the same 
for each simulation s in the simulation series: 



1 + 2<°(s) = 1 + 7if{s) + n°"(s) = Afc for each s 



(27) 



(s) 

The mean cycle length Tcyc of the path-simulations 

and transmission coefficient calculations are pro- 
fs) 

portional to the average pathlength T^J^^^ 
cepted paths: 



of the ac- 



1, 



for MD or MC, 



(s) 

^'^path path sampling and transmission calc. 

(28) 



, Tcyc are given in App. ^ and ^ and 
Tpath is the average number of MD steps of the accepted 
trajectories. Hence, we neglect the fact that Afc{s) and 
^(s) can differ for each simulation s in a simulation sc- 
ries. In fact, the equations that we will derive for TIS 
are true even if a softer assumption is valid, i.e. that 
Afc X ^ is constant for each s. In general, ^ is smaller 
than 1 as rejected pathways are usually shorter than the 
accepted ones. Some rejections are even immediateiii^ 
and do not require any force calculations. The trans- 
mission coefficient calculation has ^ = 1 as each point 
on the TS, obtained from the first free energy calcula- 
tion, with randomized Maxwellian distributed velocities 
is automatically accepted. Still, the pathlengths Tpath 
can differ. 



A. RF methods: The free energy calculation for 
the 6 — case 

As explained in Sec.m the RF methods consist of two 
independent types of calculations: the free energy and 
the transmission coefficient calculation. Moreover, there 
exist several techniques to determine these two quanti- 
ties. As the TST approximation H1U|) is exact for the 
6 = case (which basically reduces the problem to one- 
dimension), the only computational cost follows from the 
free energy calculation. Contrary, when 9 is significantly 
larger than zero we can expect that the transmission coef- 
ficient calculation is the dominant computational factor 
even though the free energy calculation becomes prob- 
lematic as well (sec Sec. lVIII . Focusing on the most dom- 
inant contributions we will therefore assume = for 



the free energy and 6 > for the transmission coefficient 
calculations. 



1. US using rectangular windows 

First, to compare the enhanced efficiency of US tech- 
niques we need to know the CPU efficiency time Toff of 
straightforward MD. Examination of Eq. © reveals that 
it simply corresponds to the calculation of the ensemble 
average of a binary function as in Eq. I|B3|) with a = 
{wm)wo- Henceforth, using our assumptions 127I28|I we 

have Teff = (il^i^Wc « = McR. e+^". It 



{wm)wq J {wm)wq 

is clear that the exponentially dependence on I3H make 
straightforward MD prohibitive for any high barrier or 
low temperature system. 

To obtain the overall efficiency for US using rectan- 
gular windows as expressed in Eq. ((T)), we first need the 

efficiency time for a single window. The principal 
result of this simulation s equals {wa)wj{ws-i)ws- The 
general approach to calculate the efficiency time for a 
composite of two averages that are obtained simultane- 
ously within the same simulation is explained in Sec. IB 2l 
In the App. 10 Eqs. I|C1IC5|) . we derive that the efficiency 
time for this single window is given by 



c 



ifr>7, 



; -4)(7+.--^)(i-e-°r) _^^ ifr<7! ^^^^ 



(l_e-o7)2 



where a = 2PH/W. Here, a is a system-specific param- 
eter, Afc is an intrinsic of the simulation, and 7 and F 
have to be optimized. If we assume that M is very large, 
we can neglect the difference of the first Wq and last wm 
windows. Following Eq. (|B2ip . the overall efficiency time 

IS given by M'^tI;^ and as M w W/{2T) is not dependent 

on 7, we can minimize rj^g as function of 7 to obtain 
also the lowest overall efficiency time. This optimum is 
achieved for half infinite windows 7 — > 00 where 



'off 



(e"r _ i)AAc (30) 
for all s. As a result, the overall efficiency time equals 

\ , .... ^3^^ 



off-(|;) (e"r-l)AAc. 



Eq. lEU has a minimum for F = 1.6a"i = 0.8W/PH 
which gives 



reft = lM{(3H)^Uc 



(32) 



The efficiency time is quadratically dependent on the bar- 
rier height H and the inverse temperature /3. Compared 
to straightforward MD this is an enormous improvement. 
The optimal window boundaries imply that the fraction 
of phase points that is sampled in the right part of the 
window is given by {ws)ws = 0.20. Hence, F should be 
adjusted to have approximately one fifth of the sampled 
points in the most right up-hill part of the window. 
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2. US using Single bias window 



In appendix [HI Eqs. HC6IC7|I . we derive the efficiency 
time for an ensemble average a = {a{x)) when it is ob- 
tained via a different ensemble using a weight function 
n-. a= (arj"i)n/(f^"^)n- Then, Tes is given by: 



- 1 



l + 2n: 
2n' 



bb 



1 



'eye 



(33) 



with b = afl^^ and c = il^^ . Note that the correlation 
functions rtc (and Tcyc for path sampling) can depend 
in principle on fl as well. As Eq. (|33|1 does not change 
whenever fl is multiplied by a single factor, we apply the 
normalization (fi) = 1. Assuming that Eqs. (|27I28|I are 
also valid in the biased ensemble, we write 



B. RF methods: The transmission coefHcient 
calculation for the 9 > case 

For the reasons explained in Sec. IV Al here we will 
study the case 9 > and, in particular, wc assume that 
PHamO > 1. Substituting a[X] = x[X] in Eq. (|B1J and 
using Eqs. (|A6I4I25I27I28|I . wc can write a general for- 
mula for the CPU efficiency time Toff for the RF method: 

TcS = ^2 J^CTcyc 

27rm/3 

« ^^e2"AAcTpath(A,(a;o)'x')^(A,(,„)_A.) (38) 

where wc used a = 2(3HRy s\n6/W ^ 1. We remind you 
that ^ = 1 for the transmission algorithm as all generated 
phase points on the surface A* are accepted and used to 
generate trajectories. In IV B II and IVB 21 we will use 
formula (|38|l to calculate the efficiency of the EC, BC2 
and EPF method. 



L a J 



(34) 



If we assume that Mc has only a weak dependency on 
O, we can minimize Eq. ()34|l by taking the functional 
derivative (5rcff/(5fi = which gives (See Eqs. HC8ICIII1 ) 



I1--I 
a 



(35) 



as optimal weight function. 

Coming back to Eq. JHJ using a = w^/, the optimal 
bias follows directly from Eq. H35|l and is given by 



' JW^ ■f A--.A<A<A-. ^^.^ 
^ 1 otherwise, 



where we used Eq. (0 and (wAf)wo = dAPA(A*) <C 1. 
Substitution of Eq. (^5)1 and a = a? = wm in Eq. 
and using that {wm^I~^) w 2[d\PA{\*)f and (O^^) w 2 
gives 



Tcff 



(37) 



Hence, a scaling behavior independent of the barrier 
height or system size can be achieved. We have assumed 
here that Nc is independent of the biasing function Q, 
which is only true for certain types of MC such as those 
in which each cycle can be generated really independently 
( hence Nc = !)• In general, MC sampling is a diffusive 
type of motion and the bias (|35I36(I should also aid the 
exploration from the top to the reactant well and back. 
Therefore optimal bias function should result in a more 
or less uniform sampling distribution, which is achieved 
when = e^^^ . Taking this into account, the efficiency 
time is a bit larger r^ff oc i?a;/dA, but still independent of 
I3H. 



BC 



In the BC algorithm the pathways are propagated only 
forward in time. Say AiTose is the longest time that the 
system takes to leave the barrier when released some- 
where at the surface A*. Hence, to have a guaranteed 
plateau in the R{t) and k{t) functions wc simply have to 
integrate Tese timcstcps so that Tpath ~ Tcsc- The second 
unknown factor in Eq. (p??^ is (A^x^)<5(a(xo)-a*)- For the 
simplified dynamics we are considering Eq. Hll() can be 
reduced to 



1 if Ay > a; and v > ^J^KeJ^, 



if Aj, < A* and v > 
otherwise 



V2A£;/r 



(39) 



with V = Xx{xo) and AE{Xy) = H—V{Xl, Xy). Following 
the same lines as in Eqs. (|23l25ll gives 



A'i+^-MA,x'(A,)e-'3^(^^^«) 



(Xxixo) X >5(A,(xo)-AJ) 



(40) 



with 



X'(A,) = 



/ d.; v^e-^-^"^-' [e{v ~ yi^) + e{v + y^)' 



/ due 



(41) 



In Eqs. H4UI41|) wc made use of the mirror symmetry along 
A*. As the Gaussian integral (|41|l has a symmetry as well 
along the v = axis, we can rewrite Eq. (j41|l as follows 



X'(A,) = 
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Substitution of Eq. in Eq. reduces Eq. (EH) to 



TrToscA/c 2c 
Teff = 5 e 



TTTescA/c ipHRy sii 

(2/3Hi?ySin6i/W^)2 



(43) 



Although, the efficiency of Eq. (|47|l has been quadrat- 
ically improved compared to Eq. H43|l . the exponential 
dependence of a = 2(3HRy s'm9/W makes this method 
prohibitive as well when is significantly different from 



zero. 



For large barriers, the exponential dependence on 
mps'mORy/W makes the method already prohibitive 
for relatively small angles 9. 



2. BC2 and EPF 

In BC2 the trajectories have to be followed until both 
forward and backward trajectories are no longer on the 
barrier. In EPF, the generated point on the surface head- 
ing toward reactant state are accepted, but assigned zero 
without further integration. Points on the surface with a 
positive velocity are first integrated backward and, then, 
integrated forward in time. This gives the advantage that 
whenever the backward trajectory recrosses the surface 
A* within a short period, this trajectory's contribution is 
assigned zero as well and the forward trajectory can be 
omitted. Hence, we have r^a'th — ^^osc and t^^^^ < Tosc 

Moreover, as S-curves are absent in this system the two 



path-functionals are identical: x [-'^l = X [-^] for all 
X. A non-zero contribution of x occurs only when both 
the backward trajectory ends in the reactant state and 
the forward trajectory ends in the product state. This 
implies that the velocity v should be positive with kinetic 
energy higher than /S.E = H — ^(A* , \y). Hence, 



BC2/EPF _ 



0{v 



2AE 



)• 



(44) 



Therefore, we can write the same equations as (|40l41f) 
with 29{v~^f^) instead ofeiv-^f^)+e{v+^f^) 
in Eq. (|41|l . which gives 



2A/^e-^^^ + erfc[y^l (45) 



We can neglect the erfc-term by invoking Eq. 1D2|I and 
omitting the terms of order 0{[(3AE]^^/^). Substitution 
of this in Eq. (|in|) gives 



3t3m 

Rye-P"e°'/a ~ 3f3m' 



and, hence, Eq. yields 



EPF < 1_BC2 _ 4TcscAAc pn 



'off 



'off 



a 



(46) 



(47) 



C. TIS: the 6) = case 

As for the RE methods, the TIS methods consist also 
of two types of simulations: the initial flux through Aq 
and the crossing probability. However, contrary to the 
RE methods, in TIS we might expect that the crossing 
probability is always the computational bottleneck. As 
Aq is placed in the low potential energy region at the foot 
of the barrier, this flux is easy to compute for all values of 
6. We will henceforth concentrate on the crossing prob- 
ability for the two cases 9 = and 6 > 0. 



1. standard TIS 

Say a(^) = VAiK+i\Xs), Ao = \* -W/2, and Am = A*. 
Eor the pathlength we write Tp^jj^ sa G{Xs — Aq)^ where 
the constants G and g will be determined later on. As 
a(^) is basically an average of a binary function, that is 
1 for the successful trajectories and otherwise, we can 
invoke Eq. (HH)) and use Eqs. (|27I28|I : 



1 - a(^) 



(48) 



Va{Xs+i\Xs) is the flux through As that reaches A^+i 
divided by the total flux through Ag-'^- for trajectories 
that come from Aq. Hence, 



{Xa;{xQ)S{Xx{xo) - Xs)h^ 



.s+1,0/ 



{X^S{Xx{xo) - Xs)hl ^ 



(49) 



where h'^^j^ equals 1 only if the backward (forward) tra- 
jectory crosses i before jii. Otherwise it is 0. Eor the 
case As < As+i < A*, we can write /ig s = ^(A(a;o)) and 
= e{X{xo) - .flKEjm) with AE = 2H{X,+i - 
Xs)/W the difference in potential energy of the two sur- 
faces. Hence, 



Pa(As+i|As) 



rOC 



2A_E/i 



- dv 'I 



-PAE 



ve 



(50) 



We take an equidistant interface separation such that 
As+i — As = = AA and As — Aq = sAA. Moreover, 



'Y' = ^a(Am|Ao) = e"^^ or, equivalently, M 



l3H/\ Ina^'*)]. This allows to rewrite Eq. lO as 

1 - 



|lna(^'r(^)'sWcCG. (51) 
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Since we know the efficiency times for each simulation 
s, we can calculate the total efficiency time of the whole 
simulation series. It is, however, important to note that 
T^g is not a constant as in Eq. (|29I30|I . but depends on 
s. This raises an additional question of how one should 
divide a given total simulation time among the differ- 
ent simulations. A logical choice would be to use the 
same simulation time for each s or to adjust the sim- 
ulation times to obtain the same relative error in each 
part. We denote the total efhciency times using these 
two strategies r'^g and t"^^ . Surprisingly, for this case the 
two approaches are equally efficient and given by 



1 



As) 



'eff 



f3H \^^W\9 Afc^G 



V|lna(-)|/ 



V 2 



(52) 



where we used Eqs. (|B19IB20|I and J^t ~ M9+^/{g^ 



1). As we show in App. IB 31 these two approaches do 
not yield the optimum efhciency. This would be attained 
when we assign each part s a simulation time propor- 



tional to 



which yields 



yielding 



1.40 
'At 



(56) 



Due to a decreasing average pathlength, TIS seems to 
have a slightly better scaling as function of H than US 
using rectangular windows H32|) . Assuming a lower prcf- 
actor for US, this will imply that the TIS efficiency ap- 
proaches the US efficiency Eq. H32() at increasing bar- 
rier heights. However, Eqs. H55I56|I break down for very 
high barriers as the average TIS pathlength Tpath can, of 
course, in practice never decrease below one MD step. 



2. TIS biasing on An 



We can exactly follow the same lines as Eqs. (I33l37f) 
which gives for the ideal biasing function 



max) ^ 2 



Va{Xb\Xa) 



if A 

max 

otherwise 



(57) 



Tcff 



(") V|lna('')| 



f3H \^ fW\9 4J^c^G 

Jg + W' 



(53) 



Minimizing Eqs. H53|l with respect to a'''-' shows that Toff 
reaches a minimum for a^"^ « 0.2. Hence, the TIS pro- 
cedure is optimized when approximately one out of five 
trajectories are successful. See the correspondence with 
US sampling lV A II where one fifth was also the optimum 
for fraction of sampling point in the right uphill part of 
the window. Using a'"^ = 0.2 in Eq. (|53|) gives 



Toft 



6.18(.^ 



W\9 {(iHf 



c 



(54) 



and 



■^^^^TjyTcff . In practice, we have found a 
1) of the pathlength on a steep 



:ff - 'cff - 4(g+l 

linear dependence {g 
barrier— and quadratically dependence (17 = 2) on a dif- 
fusive barrier—. Hence, is about 12 % and 33 % 
lower than r^g and t"^. Note that Eq. H54() has the 
same quadratically dependence on f3H as US (pl^ . One 
should realize that normally Mc^^'^'^ > Mc^'^^ and 
that A/c^^^^^ usually has a strong W dependence ex- 
cept for exceptional cases where MC cycles can be gener- 
ated really independently. Hence, the efficiency scaling of 
TIS is comparable with that of US sampling using rect- 
angular windows. As US has more flexibility to reduce 
A/c than TIS to reduce Tpath, US is probably a bit more 
efficient than TIS by a single prcfactor. 

In this specific system, the H dependence of the TIS 
efficiency is even a bit favorable than g ~ \. In the 
App. ini Eqs. l|DllD2p we derive that 



G 



2 



Wm 
H ' 



(55) 



and the overall efficiency 



Toff = ^Vc^Tpath- 



(58) 



Hence, the ideal bias function H57|) allows to obtain a scal- 
ing independent of (3H as in Eq. H37|) . Note once more 

that generally A/b^^^^^ > ^fc^^^■ As for US, if we take 
into account the diffusive behavior of the sampling, Ac 
is likely very large when the bias (|57|l is used. This bias 
favors only trajectories that reach state B, but docs not 
aid the system to climb the barrier in successive cycles. 
Therefore, a bias that generates a more uniform distribu- 
tion is more advantageous than Eq. H57|) . However, this 
does not affect the scaling dependency on f3H. 



3. other path sampling methods 

The estimation of the optimal interface separation on a 
straight slope is a bit difficult for path sampling methods 
like PPTIS, Milestoning, and FES (Note that the PPTIS 
memory-loss assumption is satisfied on the strictly in- 
creasing barrier even if the system is not diffusivoi^). An 
efficiency analysis of FFS^ using the approximation ljB8|) 
revealed that TeS is constantly decreasing as function of 
the number of interfaces. The same result would be valid 
for PPTIS. However, as correctly stated in Rcf. [sj, the 
apparent conclusion, that the computational cost can be 
made vanishingly small using an infinite set of infinitesi- 
mal spaced interfaces, is purely artifical. If one takes into 
account that there is a minimum path length (of at least 
1 MD step), also the PPTIS and FES show a quadratic 
dependence on (3H. Considering the lower path length, 
the efficiency is likely to be very close to US using rect- 
angular windows. 
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D. TIS: the e > case 

We take Aq = A* — W/ {2 cos 9) — Ryi&nO and \m = 
A* + Ryta.n9 (Sec Fig. With these outer-interface 









e / 















FIG. 4: Top view of the two-dimensional barrier and the po- 
sition of the interfaces Ao, A^/, A*, and \m- 



Xx < Ao ensures that the system is at the flat left side 
of the barrier and A^; > Xm ensures that the system has 
crossed the barrier ridge. As in Eq. (|25|l we can write 

P. (X \X)- /dA,PA(A..+r|A.;A,K/^^(^-.>^) 
i^^(A,+i|A3j „-^v(A,,A,) (5^) 



/dAj,( 



with P^(As-i-i |As; Ay) the crossing probability along the 
line Aj,. Now, the optimization of the full TIS algorithm 
would yield an extended and tedious calculation. There- 
fore, we derive an upper bound of the CPU efficiency 
time instead which is relatively easy. As in Eq. (|5(J|I . we 
can write for P^(As+i|As; Ay): 



P^iXs+i\Xs; Xy) = min[l,e 



-/3[\/(A,+i,AJ-y(A,,AJ]i 



(60) 



The additional min-tcrm compared to Eq. (|50|l is because 
the potential energy can decrease from A^. to A^+i along 
some coordinate Ay. Using equidistant interfaces with 
spacing AA = (Am — Ao)/M, we have F(As+i,Ay) — 
V{Xs,Xy) < 2AAcos6i/PF. Hence, Pr^{Xs+l\Xs■, Xy) > 



P^(A,+i|A,)>e-2A^^™«''/^ 



(61) 



We remind you that the higher a^*^ = P4(As+i|As) the 
lower Tjlg due to Eq. (|B3p . The equal sign in Eq. H61|) 
for all s, would correspond to the 9 = case with barrier 
height H' = 2MAAi7cos6'/VF = 7^(1 + ARy sin 9/ W) 
and barrier width W = W/ cos 9 - 
we can simply invoke Eq. H54() and write 



- ARy tan 9. Therefore, 



/VF'\9 (BH')'^ 
Teff <6.18(^) SP^G^N-c 



(5 + 2)2 



(62) 



which has only a quadratic dependence on PH' ^ 
(3HRySm9/W. This is far superior to exponential scal- 
ing of Eqs. H43I47|) . As a result, for high barriers and 



non-vanishing angles 9, the TIS method becomes orders 
of magnitude more efficient than the reactive flux meth- 
ods. 

Of course, one might object that the reactive flux 
methods for this 2D system improve dramatically if we 
would chose Aj^ instead of A^; as RC. The Reactive Flux 
efficiency is then again identical to the 9 = case. How- 
ever, this is exactly the crucial point. It is quite simple 
to find a proper RC in a low dimensional system, but in 
high dimensional complex systems this is certainly not 
the case. Some methods have been devised that sys- 
tematically search for RCs, but they have their limita- 
tions. The optimal hyperplane method^Si^^iS^ and the 
string method^ rely on harmonic approximations and 
on the assumptions that these hyperplanes can be de- 
scribed as a linear combination of Cartesian coordinates. 
Complex reaction mechanism, notably chemical reactions 
in solution, require a highly nonlinear function as RC. 
The inclusion of important solvent degrees of freedom is 
not an easy task. Some success has been made using 
the coordination number as RG2L2S. However, the influ- 
ence of the solvent can be more subtle. In Ref . [s^ , we 
found that electric contributions due to nearby sponta- 
neous formations of tetrahedral ordered water molecules 
can be crucial to give a last push over the potential en- 
ergy barrier. To incorporate such an effect in a one- 
dimensional RC would be an enormous task and can not 
be made without a priori insight in the mechanism. Auto- 
mated multidimensional US sampling approachesii allow 
to explore the free energy surface efficiently in a prede- 
termined set of order parameters. From the reduced free 
energy potential the most optimal one-dimensional RC 
could be estimated^. However, as the predetermined 
order-parameter space is still limited to e.g. 6 coordi- 
nates, it is still likely that important coordinates such as 
solvent degrees of freedom can be missed. 



VI. HYSTERESIS 

Up to now. we have given expressions for the efficiency 
times while treating the effective correlation Mc as a 
simple constant factor. The calculation of this factor 
is difficult as it depend on the intrinsic diffusive behav- 
ior of the MC/MD sampling itself. Fact is that Afc can 
be significant larger 1 and usually has a scaling depen- 
dence on some system parameters (like Rx,Ry, W, and 
9). Especially, the ^-dependence can be severe and will 
be discussed in this section. 

US sampling and TI constrain the system in a small 
window or on a surface that drags the system over the 
barrier. We have assumed that the time consuming step 
in the rate calculation for the 2D barrier is the trans- 
mission coefficient calculation. In fact, the calculation of 
the free energy barrier can also be very hard due to an 
improper RC. This problem is known as the hysteresis 
problem which basically results in a diverging Afc ■ Evi- 
dently, one could ask whether this problem occurs in TIS 
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as well. If this would be the case, the TIS efficiency would 
be much less advantageous as suggested by Eq. 

We will show that sampling of paths instead of phase 
points alleviates or even eliminates the hysteresis prob- 
lem entirely. The hysteresis problem does not occur in 
the 2D system we are considering, but can still persist, al- 
though still less than in free energy methods, for systems 
that have multiple reaction channels. 

Consider the potential defined by Eq. H18|) and the hy- 
persurfaces A^' , A* and Am as depicted in Fig.^ Suppose 
that we apply TI or US using small windows located at 
these surfaces. The distribution at these surfaces along 
the Ay direction is then given by 



p. 



,TI/US 



(J(Ay)J(A(x) - A,)) 
{S{X{x) - A.)) 



(63) 



where As can be cither A,,', A* and Xm- In TIS, we could 
look at the distribution of ffist crossing points with the 
As interface. This distribution is given by {v = Aa;(a::o)) 



pTIS 



{S{Xy)vS{X{x) - As)fe^.s 
(vSiX{x) - XsKJ 



(64) 





Distribution TI/US — 
Distribution TIS 


Surface: Xs' 
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\ 

\ 


Surface: ?ix ^^^^^^^^^^^y^ 




Surface: A,m / 





0.0 , 

FIG. 5: (Color online) Distributions along the Ay coordinate 
for three surfaces and the TI/US and TIS methods. 



The additional u-term in the nominator and denominator 
of Eq. (|64|l compared to Eq. H63|) is due to the fact that 
TIS looks at crossing points while the free energy distri- 
bution looks at points on the surface. A crossing point 
is a point that can cross the surface in a single timcstcp 
and because limAt-fO St^(As - A(xo))0(A(a:At) - As) = 
8(X(xq^ — Xs\X(xq^ the additional w-term appears. The 
other term g in Eq. that is missing in Eq. H63() 
ensures that not all velocities are considered. Starting 
from a; going backward in time, Aq should be hit before 
As. This implies that we can write /iq ^ = Q(v) if Aj^ < 

and /ig., 6i(u-^2(i7- y(As,Ay))/2) for A_l > 0. Sub- 
stituting this in Eq. H64() yields 



3TI/US 

As 



(A.) 



g-/3y(A,,A„) 

JdA^^WCT 



g-/3V'(A,,AJ 

/dAe-'9^'(^-^) 
(65) 



with V'(A„Ay) = V{X^,,Xy)B{-X^_) ^ HB{Xx_). Fig. 
shows the distributions of Eqs. (|65|l for three interfaces 
for the case d = 33.7, Ry = 1.5, W = 3.6 and = 9. 
At the first interface in Fig. |31 the two sampling distribu- 
tions of Eqs. H65|l are exactly the same. The distribution 
is maximized at the left side of A*. However, at the in- 
terface A* the two distributions are very different. The 
TI/US distribution has now two maxima at cither side of 
A* . As MD and most MC schemes generate a new phase 
points by making small displacements from the previous 
point, the low probability region in the middle needs to 
be crossed to sample the distribution at either side. How- 
ever, as crossing this low probability region is a rare event 
on itself, the system might not cross this region during 
whole simulation period. On the other hand, once the 
TI surface or US window is dragged over the barrier. 



the distribution is peaked at the right side of A* as can 
been seen from the shape of the distribution at Aj\/. If 
we move the surface back from Xm to A*, the system 
would be stuck again, but now at the right side of A*, 
which illustrates the hysteresis problem. Only when the 
sampling is done extensively long, both maxima at the 
surface A* can be sampled in a single simulation. How- 
ever, all measurements between two crossing events are 
correlated and basically contribute as a single uncorre- 
lated measurement. Hence, N^^^^^ becomes exceedingly 
large. 

TIS does not have this problem. The distribution at 
A* has still only one maxima. The distribution becomes 
uniform at Am- The divergence of Nq^'^ is, hence, not 
to be expected. Only when there are distinctive different 
reaction channels, ergodic sampling becomes difficult as 
for any method. Still the sampling of multiple reaction 
channels is likely more effective using path sampling than 
TI or US due to the non-locality of the shooting mova^. 
These findings and the results of IV Dl actuallv prove the 
relative insensitivity of TIS to the RC as compared to 
the RE methods. This quality has been assigned to path 
sampling methods before, but to our knowledge, this is 
the first time that this is rigorously proven for reaction 
rate calculation methods. This points out a significant 
advantage of TIS for systems for which a proper RC can 
not easily be derived. 

It is important to note that other path sampling ap- 
proaches such as PPTIS (or Milestoning) and FFS do not 
have this advantage. The PPTIS approximation fails for 
9 > Q except when the interfaces are very far apart. For 
instance, consider the interfaces As-i < As < As+i with 
As = A*. The PPTIS memory-loss assumption reveals 
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that the trajectories that start at A^-i and cross As have 
on average the same probabihty to reach Xg+i as the tra- 
jectories that start at Aq and cross A.,. This can only 
be vahd if As-i < As' (See Fig. If As-i is set closer 
to As, many trajectories that cross As-i along a lines 
Aj, < A* will not come from Ao if they are followed back- 
ward in time and the PPTIS approximation results in an 
overcstimation of the reaction rate. In turn, the large in- 
terface distances results that some of the PPTIS crossing 
probabilities are very small and can only be determined 
efficiently using TIS. 

FFS, although in principle exact, also has a problem 
for non-zero angles 9. As FFS only propagates forward in 
time using the successful paths from the previous simula- 
tion, all the simulations become correlated. This implies 
that whenever the simulation at the first interface is in er- 
ror, all other simulation results will be erroneous as well 
even if these simulations are performed infinitely long. 
This is a direct effect of the non-validity of Eq. ljB8|) for 
FFS. Coming back to Fig.jS] the flat distribution at \m 
can only be reproduced using FFS when at the previous 
interfaces (as A* and As') a significant number of path- 
ways is sampled in low-populated right tail of the path 
distribution which requires the sampling of an extremely 
large number of pathways. 



VII. CONCLUSIONS 

We have derived analytical expressions to determine 
the efficiency of different computational methods for the 
calculation of reaction rates. The efficiency has been ex- 
pressed as the computational cost to obtain an overall 
relative error of 1 when all algorithmic parameters are 
optimized. We have called this property the CPU effi- 
ciency time Teff. In App. ^ we have derived the CPU 
efficiency for very general cases. This also reveals an im- 
portant generic result of how one should divide a fixed 
total simulation time among a set of independent simula- 
tions to get the lowest overall error. After a first round of 

(s) 

simulations, reasonable estimates of t^^j can be obtained. 
Then, the minimal additional simulation time, needed in 
each simulation, to obtain the overall best performance 
can be calculated and used for a second round of simula- 
tions. This approach can be very profitable and it is not 
restricted to rate calculations. 

We have applied these rules to determine the efficiency 
of a simple 2D benchmark system that allows an analyti- 
cal approach. This offers a way to compare the efficiency 
of the different methods in a very transparent way. The 
two classes of methods that we compared arc the RF 
methods and the path sampling methods. The RF meth- 
ods consists of the calculation of the free energy barrier 
and the calculation of the transmission coefficient. Both 
for the free energy as for the transmission coefficient cal- 
culation different methods exist. For the free energy cal- 
culation we compared two approaches of US. One using 
a series of rectangular windows and one using a single 



optimal biasing potential. For the path sampling meth- 
ods, we have concentrated on the TIS technique which 
is an improvement upon the original TPS algorithm to 
calculate rates. The PPTIS approach reduces the compu- 
tational cost even more but relies on the approximation 
of 'memory loss'. As for US, we suggest that a single bias 
potential based on Amax [X] could replace the discrete set 
of interfaces. 

The 6 = case corresponds to the situation where the 
chosen RC is optimal. The TST approximation is then 
exact so that the free energy calculation is sufficient for 
the RF methods. We found an efficiency scaling equal to 
{PH)^ for US using rectangular windows. We obtained 
the same scaling rules for standard TIS. Using a single 
continuous bias potential, the US and TIS efficiency can, 
in the optimal case, become independent of the barrier 
height and temperature. However, knowing the optimal 
bias basically implies that one already knows the answer. 
US using several windows and standard TIS are more ro- 
bust approaches if no a-priori knowledge of the system 
is available. This shows that TIS and US compare very 
well in efficiency for all barrier heights and temperatures 
and that the difference can only rely in a single prefac- 
tor. It is likely that US is a bit better than TIS that 
has to be faithful to the natural dynamics of the system. 
US has more ffexibility to optimize the method such as 
choosing the best MC moves that minimizes the number 
of correlations. 

When 9 > 0, the chosen RC, A^,, is no longer opti- 
mal. In contrast to the 9 = case, the principal compu- 
tational effort of the RF method lies in the calculation 
of the dynamical correction. We showed that the BC 
method scales as ^ exp{4:/3H sm9Ry/W), while the BC2 
and EPF methods are quadratically more efficient. The 
EPF method is, however, a bit more than a factor 2 more 
faster than the BC2 method. The exponential depen- 
dence on f3H sin 9 indicates that a small deviation from 
the optimal RC can have dramatic consequences for the 
efficiency of the RF methods if they are applied to high 
barrier or low temperature systems. In contrast, the TIS 
efficiency scaling is only ~ {(3H sin 9 Ry /W)"^ . We also 
discuss the potential problem of hysteresis in US and TI 
when a non-optimal RC is chosen and why this problem 
does not occur for TIS for the 2D barrier. This gives an- 
other evidence that the TIS method is less sensitive the 
choice of RC. 

The advantage that path sampling does not require a 
RC has been advocated many times. However, although 
this statement is quite evident if the main object is to 
sample a representative set of reactive trajectories, it 
is not so evident for the calculation of reaction rates. 
The calculation of the reaction rate still requires a RC 
(the only exception we know of is the method proposed 
m Ref. hj using the pathlength as transition parame- 
ter). However, we are now the first to prove that a path 
sampling-based reaction rate calculation method, using 
TIS, is potentially orders of magnitude faster than the 
RF-based methods whenever the right RC can not be 
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determined. The reason is because TIS uses an impor- 
tance samphng technique to calculate the dynamical fac- 
tor. The importance sampling of the RF methods only 
concentrates on the free energy. Therefore, whenever 
the dynamical factor is low (e.g. due to a wrong choice 
of RC), these methods automatically run into problems. 
The main question remains whether it is more profitable 
to search for a good RC and use the RF methods, or take 
a simple order parameter (non-optimal RC) and use the 
TIS method. This question has not an easy answer and 
depends on the complexity of the system. 
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grateful to Rosalind Allen who made useful suggestions 
to improve the first version of this paper. This work 
was support by a Marie Curie Intra-European Fellow- 
ships (MEIF-CT-2003-501976) within the 6th European 
Community Framework Programme. 



APPENDIX A: GENERAL DEFINITIONS 

1. Ensembles averages in phase and path space 

We denote with x a phase space point {r, p} where r 
are the Cartesian coordinates and p the momenta. The 
examples presented here consider only 1 particle in a two- 
dimensional potential, but x,r,p are in general multidi- 
mensional vectors. The distribution of x is given by the 
probability density p{x). In case of Boltzmann statistics 
p{x) oc e^^-^'^^'^ with E{x) the total energy at phase point 
X and f3 = l/{kBT) with T the temperature and ks the 
Boltzmann constant. Suppose a is the value of a quan- 
tity we want to compute. In many cases such a quantity 
equals the expectation value of a certain observable. We 
denote this observable with a{x) which is a function of 
X. Then the expectation value or population mean (a) is 
given by 



{a{x)) 



J dx a{x)p{x) 
j dxp{x) 



(Al) 



and a = (a) for this specific case. Path sampling 
simulations require a more extensive formulation espe- 
cially when stochastic dynamics is applied. The dis- 
crete representation of a path is the most convenient, 
where a path X is defined as a set of + + 
1 successive phase points that determine the system 
at intervals of At along a certain trajectory: X = 

{a^-r''Ati X{l^T'')Ati ■■ -1 X-At; Xq; XAt, ■ ■■ , X+rfAt, }■ 

TIS, the backward r*" and forward r-^ time indices are not 
fixed, but depend on when a certain interface is crossed. 
The weight of the path is given by the initial distribution 
at time t = and the probability densities corresponding 



to the history and future of the path: 

P[X] = p{xo)Y[Pn{X(i-.,)At ^ X-iAt) 
1=1 

X n Pnix{i-l)At XiAt) (A2) 

i=l 

where Pn{x — )■ y) is the probability that the natural dy- 
namics of the system brings you from x to y given you 
are in a;, and Pnix *— y) is the probability that if you are 
in X you came from y in the past. As shown in Ref. |l4l |. 
whenever the system is in a steady state, Eq. (|A2ll is 
identical to 

P[^] = PiX-rOAt) Yi PniX{i^ri>-l)At X(^,^r'')At) 
i=l 

(A3) 

which corresponds to the path weight as expressed in the 
original TPS papers^ with the only difference that the 
starting index is — r** instead of 0. Now, if a[X] is a 
function defined in path space, the population mean is 
given by 



i[X]) = 



JVXa[X]P[X] 
JVXP[X] ' 



(A4) 



with VX = nr=_rf dx^At and P[X] given by Ea. l|J2|) . 

In the following, we use g as a point that is defined in 
either phase or path space, i.e. q can either represent x or 
X. Besides the simple ensemble averages (HU and (HU, 
we can also define the biased ensemble average (a) q using 
a weight function ^(q). This biased ensemble is defined 
as 



(a(g)»(g)) 



(A5) 



In practice, sampling the biased ensemble {■■■)n 
by means of MD simply requires the addition of a 
term — (Infi)//? to the potential of the system. In 
MC, this is achieved by a change of the acceptance- 
rejection criterion from min[l, p(a;^"-')//9(a;^°^)] to 
min[l,p(x("))17(a;("))/(p(a;(°))f}(j;(°)))] with a;(") and 
a;(°) the new and old generated MC points^. 

Ensemble averages in path space (a [A]) can depend on 
the type of dynamics (hence, on the hopping probabili- 
ties pn). To specify this dependency, explicitly, we use 
(. . .);p„, which allows to generalize these path ensemble 
averages to arbitrary dynamics. The hopping probabil- 
ities of a simulation method Ps{x — s- y) can be of any 
type, for instance Monte Carlo, and do not need to be 
related to the natural dynamics of the system. There- 
fore, these ensemble averages are annotated by {. . .)._p^. 
When ;p is not specified, we assume that the natural 
dynamics is applied or that the property is independent 
of the hopping probability p. Hence, in our notation 
{a[X]) = (a[A])i = {a[X]),,„ - {a[X]),,,^. 
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2. Standard deviation, variance, covariance, and 
correlation 



The standard deviation of d is defined as 



aa = («))') = 

= Viia-a^) = V(a'> - a'- (A6) 
Related to the standard deviation is the variance of d 

Var(a)^((a-(a)f)=a2 (A7) 

and the covariance of two functions d and b 

Cov(a,6)EE((a-(a))(5-(6))) 

= {db) - ab (A8) 

with Cov(a, a) — Var(a). A simulation s, which can ei- 
ther be MC, MD, or TPS/TIS generates a set of points 
{^Oj 9ij 92, ■ • ■ , 9n} in either phase or path space. Each 
point qi is simulated with a certain probability Psiqi) and 
the chance that after qi another point qi^i is sampled is 
given by Psiqi — > Qi+i)- We denote Ui = a((ji). Now, the 
sample mean d(n) for a simulation of length n is given 
by 



a{n) 



1 " 

77 ^ 



(A9) 



Eq. IIA9|l converges to the exact value, a{n) = a, in the 



limit 



oo. The standard deviation in the mean a 



a(n) 



is defined as the standard deviation in the set of points 
{a(i)(n),a(2)(n),a(3)(n),a('')(n), . . .} that is obtained af- 
ter performing a large number of independent simula- 
tions, s = 1,2,3,..., each of length n and with result 



(n). Hence, we can write 



^ n n 



n-l 



.1 J=l 



2—1 



(AlO) 



Using time translation invariance, we can show that 
((flj - af).p^ = ((ao - a)2).p^ = ((a - af) and {{a.i - 
a){aj — a));p^ = ((oq — a){aj-i — a))-p^. Now, we assume 
an exponentially decay in correlation for large I = j — i 
which yields for large n: 



2 



^ oo 

-{((a - af) + 2 ^((ao - - a)),p^ } 



which gives 



2E 



((ao - a){ai - a)).p^ 



1=1 



{{d-ay 



} (All) 



(A12) 



In Eq. l|IT2|) . is the correlation number or the inte- 
grated auto-correlation function, which is a special case 



of n°^, defined as 



with 



(A13) 



1=1 



C^,^1)^(S^1^^}(^LZ3^, (A14) 
{{d-a){b-b)) 

If two functions a and b are uncorrelated, (a6) = (a) (6) = 
a6. Hence, if the simulation data are not correlated 

Psiqi qt+i) = Ps{qt+i) and ((a; - a)(aj+; - a)) = 
{{a, - a)){{a,+i - a)) = and = 0. 



3. Statistical errors and propagation rules 

Let Aa be the absolute error of a quantity a and let 
6a = Aa/a be the relative error. For a sequence of mea- 
surements {ao, ai, . . . , a„}, the absolute error in the av- 
erage d(n) is defined as the standard deviation in the 
mean: 

Aa(n) = cra(„). (A15) 
Using Eq. (|A12() . the absolute and relative errors are 



Aa(n) = da 



1 + 2n2 



Sa(n) 



aa_ ll + 2nT 
a 



n 



(A16) 



This shows that it is of eminent importance for the effi- 
ciency of the method that the trajectories generated by 
the computer algorithm minimize the correlation num- 
ber Tic as much as possible. In MC, this gives rise to 
conflicting strategies. In general MC algorithms select 
a new phase point by making a random displacement 
from a previous point, which is then accepted or rejected. 
When the displacement is small, the new point resembles 
strongly to the old one which can be a source of correla- 
tions. On the other hand, if the randomized displacement 
is too large, it is likely rejected after which the old phase 
point is counted again. Therefore, a high rejection prob- 
ability also introduces correlations. The optimum max- 
imum displacement is a compromise between these two 
effects. The creation of a single step in path sampling 
is much more expensive than standard MC or MD, but, 
on the other hand, the correlation number Uc is usually 
much lower. 

Suppose a quantity a is not equal to the expectation 
value of a single operator, but, for instance, depends on 
two other quantities b and c via a function f: a ~ f{b, c). 
The error in a can then be determined using the error 
propagation rules. Say b(rn) and c(n) are the approxima- 
tions of b and c after m and n simulation cycles respec- 
tively. The resulting error Aa(m, n) can then be derived 
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as follows. As (b{m) — b) and (c(n) — c) arc small for 
large values of m and n, we can invoke following Taylor 
expansion: 

a{m,n) = f{b(m),c(n)) 

^ f(b,c) + ^{b{m)^b) + ^{c{n)~c). 

(A17) 

As in Eq. (|A15p the error Aa(m,n) equals the standard 
deviation of the mean a(m, n). Hence, 



/\a'{m,n)^{{a{m,n)^af),,^ = {^^)\{b{m) - bf) 



dcj 



(A18) 



Substitution of Eq. IjASp . (|A15|) and l|A12|l yields 

+ 2^^CoY {b{m),c{n)). (A19) 

Then, from Aa = /(&, c)(5a, A6 = 6(56, Ac = cSc we get 
the propagation rule for the relative error 

6a^im,n) = [b^ / f]' 6b^ (m) + [c^/f]'ScHn) 



dc 



+ ^[%%/f]Cov(b{m),c{n)). (A20) 

Eqs. l|IT9|) and (|X20)l arc the most general rules for the 
propagation of errors when a depends on two quantities 
6 and c. It can straightforwardly be generalized to more 
arguments such that a = f{b, c,d, . . .). 



APPENDIX B: CPU EFFICIENCY TIMES 

1. CPU efficiency time for a single ensemble 
average 

Let us define Tcyc as the average CPU time to perform 
a single MC, MD or path sampling cycle. Moreover, we 
denote Tgim = riTcyc as the total CPU simulation time 
after n cycles. We introduce now the CPU efficiency time 
Teff(a) to obtain a relative error Sa equal to one. Putting 
i5a = 1 in Eq. ljA16|) and using n = Tsim/Tcyc = Tcs/Tcyc 
gives 



Tcff(a) 



^ (l + 2nrK: 
a / 



(Bl) 



This is the general CPU efficiency time to obtain a rela- 
tive error of 1 in a single simulation average. This CPU 
efficiency time is a characteristic property of the quan- 
tity a (actually of the whole function a) and the simula- 
tion method (via the ps hopping probabilities and Tcyc). 



Whenever Tcff(a) is known the absolute and relative er- 
rors after a simulation period Tsmi are directly given by: 



Aa = a^roff(a)/Tsi,n, 6a = \/T\^a{ayT~. (B2) 

If a is a binary operator such that a{q) is cither 1 or 0, 
then a? = a in Eq. (|A6p and Eq. IjBip equals 

rcff(a) = ^(l + 2<")rcyc. (B3) 



Therefore, the Tcff (a) becomes very large for small values 
of a and, hence, an accurate evaluation becomes prob- 
lematic. 



2. CPU efficiency time for a composite quantity 

The CPU efficiency time for a composite value a = 
f{b,c), where 6 = (6) and c = (c) are determined simul- 
taneously in a single simulation, can be derived following 
the the same lines as in Eq. HA11|1 for the covariance 

Cov(6(n),c») = -{ ^((6,; - 6)(c, - c));^^ 

n n 

+ E E - b){c, - c));p, + ((c, - c){b, - 6));pj} 

i=l 

«i{((6-6)(c-c)) 



1+ 



E 



((6o - 6)(q - c)).p^ + ((co - c){bi - b)),p^ 



1=1 



((6-6)(c-c)) 



} 



-Cov(6,c) 
n 



(B4) 



Here, we inserted Eq. (|A8p and (|A14p in the last line. 
Substitution of Eas. l|B2|) and ||B4|) into Eq. HA20|I and 

using n = Tgim/Tcyc gives 



^a=[6-//] 
-df df 



\ "J IfV 



+ 2[S;g://']Cov(&, c){l + nl^ + }-^. (B5) 
Taking Sa^ ~ 1 and Tgim ~ Tcff (a) directly results in 



reff(a) ^ [6^//] Vcff(&) + [c^//] Veff(c) (B6) 



where Toff (6) and Tcs{c) are given by Eq. (|B1|I . For ex- 
ample, if a = b''cJ , the efhciency time of a equals 



Tcff (a) = i^rcff(6) -I- j^rcff(c) 



2^Cov(6,c){l + n!;^ + < 



J ' eye • 



(B7) 
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3. CPU efficiency time for a series of simulations 

Suppose a — f{b,c) and b{m) and c{n) are obtained 
via two different simulations. Hence, n and m are not 
necessarily the same. In the following we assume that the 
different simulations are uncorrelatcd. Thus, we assume 
that for two different simulations (1) and (2) the following 
holds 



Cov{b{m),c{n)) = {{b{n) 



&)(cM -c));p^,, ,) = 0. 

(B8) 

Here, the subscript ;Ps(i,2) indicates that the ensemble 
average can depend on how the two simulations are con- 
nected. Assumption ljB8|l is, for instance, true for US 
when two independent simulations are performed using 
two overlapping windows. It also holds for TIS where 
the outcome of an interface sampling simulation at a cer- 
tain interface is independent of the result of the previous 
interface simulation results. Eq. (|B8|I does not hold for 
the most common implementation of the reactive flux 
method. In this approach, the importance sampling to 
determine the free energy barrier is simultaneously used 
to obtain a representative set of configuration points at 
the TS dividing surface. These configurations with ran- 
domized Gaussian distributed velocities initiate the dy- 
namical trajectories that determine the transmission co- 
efficient^i^. Moreover, at variance with TIS, Eq. (|B8p is 
not true for the Forward Flux Sampling (FFS) method, 
that was devised by Allen et al^^^. Here, the result of 
the interface sampling at one interface depends on the 
results of all the previous interfaces. The importance of 
Eq. HB8|) is further discussed in Sec. IVII 

If the total simulation time Tsi,„(a) = Tgimib) + 
Tsimlc) = niTcydb) + nTcyc(c) is fixed, we still have some 
freedom in choosing n and m or, equivalently, choosing 
Tsim(fe) and Tsini(c). First, by substitution of Eqs. HB2p 
and HB8|) into Eq. HA20|I we obtain 



Tcsib) 



Toff(c) 



(B9) 



1 

(&) ' '- dc' ' rsi,n(c) 

A logical approach, although not the optimum, is to give 
each simulation the same simulation time. We will denote 
the efficiency time that results from this strategy r'^g^a). 

Taking Ts\ni{b) = Tsim(c) = ^Tsun{a) for T^ff(a) = Tsim{a) 

and Sa^ = 1 gives 



(BIO) 



Alternatively, we could try to obtain the same error in 
each simulation. The corresponding efficiency time will 
be annotated as r"ff. Then, we need to use simulation 
times proportional to cx r^ff (&), Teff (c). Hence, we take 

Tsim(^^) = Tsiin(a)Tcfr(6)/[Tcfr(6) + Toff(c)] and Tsim(c) = 

Tsim(a)Tcff(c)/[Tcff(&) 4- Tcff{c)] in Eq. (|S^ which results 
in 



f(a) = (r^sib) + Tcff(c)^ 



(BID 



In order to determine the lowest TeS, we add Lagrange- 
multipliers constrains to Eq. (|B9|I in order to fix the total 
simulation time Tsim(a) 



nib),Ts^mic),r]] = [b-g^/f] 
r df , „i2 T^ff(c) 2/ / X 



.(c)} 
(B12) 



and minimize Eq. (|B12() with respect to all its arguments. 
Taking the derivative to Tsini{b) gives 



Therefore, 



.ib) 



Tsim{b)y 



(B13) 



(B14) 



We have put the absolute signs | • | as the simulation times 
needs to be positive. The same relation holds for rsiin(c). 



.(c) 



dc 1] 



(B15) 



We can sum up Eqs. ljB14|l and (|B15|) and use Tsi„i(a) = 
'''sim(^) + ''"sim(c) which givcs thc solution for rj 



1] 



(B16) 



Then, substitution of Eq. (|B16() in Eqs. I|B14IB15|I com- 
pletes the equations for Tsh-n{b) and Tsi„i(c). Substitution 
of these two equations in Eq. ljB9|) results in 



off («) = (I [b^/f] I + 1 [c^/f] I v^y. 

(B17) 



df 



The efficiency time Tcff (a) of Eq (|B17|) is always strictly 
less than or equal to T^ff (a) and r^jj (a) of Eqs. (|B10IB11|I . 
These CPU efficiency times are straightforwardly gen- 
eralized to simulation series of any number. Suppose 
that the final desired value a is obtained by a = 
f{a^^\ a^'^K . . . , a'^^-*), where a'-*' refers to the exact value 
that should be produced by simulation s and M is thc to- 
tal number of independent simulations that are needed to 
determine a. Then, the CPU efficiency times Eqs (|B10t 
IbT7|) yield 



•^ofr(a) 



M 



M M 



s=l 
M 



Tcff(a) = (XI [' 



s=l 



df 



/f\ 



If] 



(B18) 
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where r^g = reff(a('*^) is the efficiency time of simulation 
s. For example, in many methods, the final value a is 



given by a product of simulation results: a = YitLi ^ 
Then, | [a^^^ //] I = 1 for any s after which Eqs. (|BT8|) 
become 



M 



and 



(B19) 



(B20) 



The same Eqs. (|B19IB20(I are valid for the case a = 
ULi a^'^/ ntLj+i with any J e [0, M]. 

Whenever t^^ is the same for all s, Eq. (|B19|1 and 
(IB20p become identical and equal to 



(B21) 



APPENDIX C: EFFICIENCY OF US 
TECHNIQUES 

(s) 

We will derive the efficiency time T^f^ of a sin- 
gle window in US as depicted in Fig. The factor 
that needs to be computed in simulation s is a^'*-' = 
{ws)wj{ws-i)w, = b/c. Here, b = b = Ws, 

c = {c)wsi Eind c = Ws-i- Then, we can use Eq. ljB7|) 
with z = 1 and i = — 1 



■s' = Tesib) + Teff(c) - ^Cov(6, c)(l + + )re, 



(CI) 



and applying the assumptions Eqs. H27I28|) gives 

T^J = Teff (6) + Teff (C) - ^Cov(6, c)Mc ■ (C2) 

Via Eq. (HHJ, Eq. and Eqs. |77l^ we get 



Toff (6) = ^—r—^C, Toft(c) = -Mc, 

c 

Cov(fe, c) — {■Ws-iWs)w!, — be. 



(C3) 



Then using Eq. Q and ((TH|l and writing A^, = A* - 
and a = /32iJ/VF we arrive at 



/.f^'dAe-^ 1- 




/(f_+^,dAe-- e-- 




_J(s-i)r 




J(,s-i)r '^-^^ 


gar _ g-a7 



(C4) 



Moreover, (u's-iWs)iVs is only nonzero in case 7 > F. 
Hence, 



0(7 - F) 



(s-i)r+7 

sT 



dXe 



sr+7 
(s-i) 

1 _ e-"(''~r) 



rsi+7 i\ 



g--Q7 



(C5) 



Substitution of Eqs (|C4IC5|) in Eqs. (|C3p and, after that, 
substitution of Eqs. (jHH in Eq. yields Eq. (|^ . 

In order to derive the efhciency time for an ensemble 
average a = {a(x)) when it is obtained by a weighted 
ensemble via a = {dn~^)n / {fl~^)n, we write again a = 
b/c with b = (6)fi, b ~ afi^-^, c = (c)o, and c = fi^-'^. 
Applying Eq. ljA5|l gives b = and c 



_ - - (ny- Then 
applying Eqs. (|A6() . (|A8p using the ensemble (. . .)o gives: 



Cov(6,c) = (&c)o-(6)n(c)o 



Substitution of Eqs. (|n6)l in Eq. llBlT) yields 



(C6) 



Toff(&) 



- 1 



reff(c) = -1) 



1 + 2n: 
1 + 2n 



bb 



'eye, 



(C7) 



Substitution of Eqs. ^C6IC7p in Eq. (ETJ yields Eq. 

To obtain the optimal biasing function, we add a La- 
grange multiplier to Eq. ^ to fulfill the normalization 
constraint 



Toff (a) = I 



(an-') 



a) ) 



-rf[{n)-l]. 



(C8) 



The functional derivative of {af{fl)) to for any opera- 
tor a and function / is given by 



S{af{n)) _ a{x)f(9)e-^^^^ 



sn 



-l3E{x) 



a{x)f'{n)P{x) (C9) 



with P{x) = exp{—PE{x))/fdxcxp{—PE{x)). Hence 
taking Stcs/SH = in Eq. ljC8p results in 



1 a' 



a^n^{x) 02 (x) ail^ix) 



N'c+Tl''}P{x) = 
(CIO) 



which has as solution 



n'{x) = ^ 



a{x) 



\ a 



(CU) 
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rf' can be obtained, if necessary, via the normalization 
requirement {Q) = 1. Moreover, from Eq. (|C11|I follows 
directly Eq. H35|l as should also obey il{x) > for 
each X. 



with At the MD timcstep that is required to express Tpath 
as an integer representing the average number of discrete 
timesteps. In Eq. (|D1|I . we have first used the approxi- 
mation^ 



APPENDIX D: TIS PATHLENGTH 

The TIS pathlength H55|l can be obtained as follows. 
Say V is the velocity at the foot of the barrier (Aq) at a 
time t — 0. The 'returning time' is obtained by solving 
following equation X{xt) = Xq + vt + -^t^ = Xq + vt — 
j^t^ Ao, which has a solution for = L{v). 

Within the path ensemble s all trajectories should cross 
Ag. Hence, at the foot of the barrier the kinetic energy 
imu^ must be larger than E = 2H{Xs — Xo)/W. There- 
fore, for the average pathlength we can write 



■fc(x) « ^ " (D2) 

VTT X+ yjx'^ + 4/7r 



_ Wm 2y/PE + ^e^^ erfc(^/^) 



Wm ^ y \J P(4:+pE7.)J ^ w and, then, neglected the 0{E-^/'^) terms. Using E = 

AtH AtH ' 2i?(As - Ao)/VK and 

'''path — Gi^Xg — Ag)^ results in 

(Di) Eq. inni). 
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APPENDIX E: LIST OF SYMBOLS 



T temperature 

/? = 1 /T inverse temperature 

ks Boltzmann constant 

r configuration point 

p momentum point 

X = (r, p) phase point 

m particle mass 

V{r) potential energy of r 

E{x) ~ V{r) +p^/2m total energy of x 

p{x) equilibrium distribution, p{x) = e"'^^^^) for Boltzmann statistics 

Xt phase point at time t 

At MD timestep 

X path consisting of discrete timeslices: {aJ-r'Atj ■ • ■ ! ^o] ■ ■ ■ i ^^t/a*} 

[X] , [X] the start and end time index of path X 

Pnix — > y) the probability density to go to y from x in one timestep by MD 

Pnix <— y) the chance that when you are in a;, you came from y one timestep before 
Ps{x — > y),Ps(x <— y) same hopping rates for two consecutive simulation cycles using MD/MC 

P[X] weight of the path X 

q point in either phase or path space 

Qi phase/path point generated after the i-th simulation cycle 

(...) ensemble average 

(. . .)o weighted ensemble average using weight function Cl(x) 

a exact value of an observable 

a(q) generic operator that determines a 

Ui = a((?j) function value of the i-th simulation cycle 

a(n) average function value over n simulation cycles 

Aa(n) absolute error in a after n cycles 

Sa(n) relative error in a after n cycles 

(7a standard deviation of the distribution {a^} 

(7a(„) standard deviation of the distribution {a{n)} 

Var(a) variance of d 

Cov(a, b) covariance of d and b 

Ca,b{l) correlation function 

n"'*' correlation number 

Afc effective total correlation 

s index of a simulation in a simulation series 

a'-*^ exact value of an observable obtained from simulation s 

Tcyc average duration of a simulation cycle 

Tsim total simulation time 

Tpath average path length in a path sampling simulation 

Tosc maximum time needed to leave the barrier region 

Teff (a) lowest computational cost needed to determine a with a relative error equal to one 

^ average ratio Tcyc/rpath 

k reaction rate 

TZ unnormalized transmission coefficient 

k, TZ time-dependent rate and transmission functions 

x[X] recrossing correction functional 

A(a;) reaction coordinate 

P(A') probability density to be on the surface {a;|A(x) = A'} 

Pa(A') probability density to be on the surface {a;|A(a:;) = A'} given you are in A 
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F{X) = — ln(P(A))//3 free energy profile along A 

A* transition state value or the maximum in F{X) 

Xs value defining interface s: {a;|A(x) = A^} 

Xa = Ao interface defining stable state A 

As = Aj\/ interface defining stable state B 

M total number of simulations used in a simulation series 

7'a(A|A') crossing probability from interface A' to A 

RxjRy dimensions of the reactant well 

W width of the barrier 

H height of the barrier 

Xx assumed reaction coordinate in the 2D system 

Xy important other degrees of freedom in the 2D system 

A_L unknown ideal reaction coordinate 

9 angle between A^: and Aj^ giving the deviation from the optimal RC 

r,7 dimensions of rectangular windows used in US 

Ws{x), Ws{x) block functions defined by Eq. Q 

pTi/us^^^-j sampling distribution along Ay at the surface As when using TI or US 

P^^^{Xy) sampling distribution along Ay of first crossing points with surface As when using TIS 

^ (s) 

g, G exponent and pre-exponential factor the assumed behavior of t^J^^yi 

hj\^{x) history dependent function that measures whether x was more recently in A than in B 

hJl Ax) history dependent function that measures whether x has crossed A^ more recent than A, 

hij{x) future dependent function that measures whether x will cross Xi before Aj 

(j)A flux function through Xa equal to S{X{x) — Xa)X9{X) 

I 



APPENDIX F: LIST OF ABBREVIATIONS 

BC Bennett-Chandler formalism 

BC2 history dependent BC 

EPF the effective positive flux 

FFS forward flux sampling 

MD molecular dynamics 

MC Monte Carlo 

RF reactive flux method 

RC reaction coordinate 

TI thermodynamic integration 

TIS transition interface sampling 

TPS transition path sampling 

TS transition state 

TST transition state theory 

US umbrella sampling 
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